*			  ==========================================
*						     Stata公开课 
*  
*                           软件及计量基础   
*             ==========================================

*                    	 主讲人：候丹丹
*                   	 主办方：连享会（www.lianxh.cn）
*           课程主页：https://gitee.com/lianxh/stataopen
*                     :: 课件下载，答疑等 ::

*			       ============================
*					 	第四讲 回归分析（一）   
*                  ============================

**-注意：执行后续命令之前，请先执行如下几条命令

  global Path 	"`c(sysdir_stata)'ado\personal\open4_regress"  //定义课程路径 
  global D		"$Path\data"		//范例数据
  global Out	"$Path\out"		//范例数据  
  cd "$D"     
  set scheme s2color  //设定图形格式为默认格式
  
  *-note: 
  /*	`c(sysdir_stata)'是一个暂元，里面存放了Stata的安装路径：
	输入 sysdir 后显示的第一个文件路径。
	例如, 我的 stata17存放于D盘study文件下, 
	所以, `c(sysdir_stata)' = D:\study\Stata17\     	*/

	

  * ===== 学习过程遇到问题怎么办？====
  *
  * ❓ 提问前，先用 -lianxh- 或 -songbl- 命令搜索相关资料和推文
      ssc install lianxh,replace
	  ssc install songbl,replace
      lianxh DID 倍分法
      songbl 面板
	  help lianxh 
      help songbl	  
  * ❓ 亦可，前往课程主页 https://gitee.com/arlionn/PX --> wiki 
  *   查看 FAQs，通常都能找到答案。 
  *
  
*----------------
*   本讲目录  
*----------------
* 4.1 基本统计量
* 4.2 相关系数
* 4.3 OLS估计原理和估计方法
* 4.4 结果的呈现和输出





 
*			  ==========================================
*						     Stata公开课 
*  
*                           软件及计量基础   
*             ==========================================

*                    	 主讲人：候丹丹
*                   	 主办方：连享会（www.lianxh.cn）
*           课程主页：https://gitee.com/lianxh/stataopen
*                     :: 课件下载，答疑等 ::

*			       ============================
*					 	第四讲 回归分析（一）   
*                  ============================
*					    -4.1- 基本统计量
*  
*-4.1.1- 基本统计量
	
	shellout "sum1.pdf"
	
	*-	sum
	sysuse auto, clear
	summarize
	sum price mpg weight
	sum price mpg weight, detail
		
	bysort foreign: sum price mpg weight  //分组汇总
	
	*-	tabstat
	h tabstat
	global xx "price weight mpg rep78"  
	
	tabstat $xx, ///
		stat(N mean sd min p50 max) col(s) format(%6.3f)
		
	tabstat $xx, by(foreign)   ///
    stat(N mean sd min p50 max) nototal long  ///
    col(stat) format(%6.1f)          //分组统计
	
*-4.1.2- 组间均值差异

    shellout "sum2.pdf"        
	ttable2 $xx, by(foreign) f(%6.2f)
	ttable3 $xx, by(foreign) f(%6.2f) tvalue	//组间均值差异附带t值
	ttable3 $xx, by(foreign) media  //组间中位数差异	
	
*-4.1.3- 结果输出

	sysuse auto, clear
	
	*-	输出描述性统计分析
	logout, save("$Out\Table1_sum01") excel replace:  ///
		tabstat $xx, stat(N mean sd min p50 max)  ///
		col(s) format(%10.3f)
			
	logout, save("$Out\Table1_sum01") excel replace:  ///		
		tabstat $xx, by(foreign)   ///
		stat(mean sd min p50 max) nototal long  ///
		col(stat) format(%6.1f)   
	
	sum2docx $xx using "$Out\sum2docx.docx", ///
		replace stats(N mean(%9.2f) sd(%9.2f) p25 median p75 min max)    ///
		title("表1 描述性统计表")  ///
		note("this is auto's Summary Statistics")
		
	*-	输出组间均值差异检验	
	putdocx clear
	t2docx $xx using t2docx.docx,   ///
        replace by(foreign)   title("表3 t检验")	

	

	









*			  ==========================================
*						     Stata公开课 
*  
*                           软件及计量基础   
*             ==========================================

*                    	 主讲人：候丹丹
*                   	 主办方：连享会（www.lianxh.cn）
*           课程主页：https://gitee.com/lianxh/stataopen
*                     :: 课件下载，答疑等 ::

*			       ============================
*					 	第四讲 回归分析（一）   
*                  ============================
*					     -4.2- 相关系数
*  
*-4.2.1- Pearson 相关系数  

	sysuse nlsw88,clear
	global x "age grade wage hours ttl_exp tenure"
	corr $x
    pwcorr $x      //缺陷: (1)小数点后两位为宜; (2)没有标注显著水平;
    pwcorr $x, star(0.05) //小数点后两位不易调整;
  
  *-自编命令
    pwcorr_a $x
    pwcorr_c $x, star(0.05) format(%7.2f) //比较符合多数期刊的要求

*-4.2.2- Spearman 相关系数

	spearman $x, star(0.05)

	shellout "Spearman相关系数.pdf"

	/*  注: 
	(1) Pearson 相关系数：
			① 反映两变量之间的线性关系；
			② 要求两变量服从正态分布；
			③ 要求两变量均为连续变量，对极端值敏感。
			
	    spearman 相关系数（非参数检验）:
			① 反映两变量之间的单调关系；
			② 适用于连续和离散序数变量，对极端值不敏感.
		
		连续变量, 正态分布, 线性关系。二者均可, Pearson 相关系数较好;
		上述任一条件不满足，用 spearman 相关系数，不能用 Pearson 相关系数。   */
		
  
*-4.2.3- Spearman 和 Pearson 相关系数矩阵的合并呈现
  
	corsp $x, format(%7.3f)
	* Pearson 相关系数位于下三角;Spearman相关系数位于上三角		
		

*-4.2.4- 输出结果
	logout, save("$Out\Table2_corr") excel replace: ///
            pwcorr_a $x
			
	corr2docx $x using corr2docx.docx, ///
		replace star(* 0.05) fmt(%4.2f) ///
		title("表2 相关系数矩阵")		//Pearson 相关系数下三角;Spearman相关系数上三角
			







*			  ==========================================
*						     Stata公开课 
*  
*                           软件及计量基础   
*             ==========================================

*                    	 主讲人：候丹丹
*                   	 主办方：连享会（www.lianxh.cn）
*           课程主页：https://gitee.com/lianxh/stataopen
*                     :: 课件下载，答疑等 ::

*			       ============================
*					 	第四讲 回归分析（一）   
*                  ============================
*					-4.3- OLS估计原理和估计方法
*  
*-4.3.1- OLS估计原理
	
	*-变量之间的关系
  
	do "R1_reg_Fig0a.do" 
  
	*-OLS 的游戏规则: 残差平方和最小
  
	do "R1_reg_Fig1.do" 

*-4.3.2- OLS估计假设

    * H1: 线性假定

    * H2: X 是满秩的，i.e. rank(X) = k
  
    * H3: 干扰项的条件期望为零（模型的设定是正确的）  
    *    E[Xu] = 0  or   E[u|X] = 0 	

*-4.3.3- OLS估计及结果解读

	sysuse "auto.dta", clear
	regress price weight 

	*-拟合值和残差   || 详见 R6_predict.do
	
	predict price_fit, xb  // 拟合值, xb 选项可以省略,默认
	gen price_fit2 = _b[_cons] + _b[weight]*weight //手动计算
	
	predict e, residual    // 残差, residual 选项是必须的, 可以简写为 r
    gen e2 = price - price_fit //手动计算
	
	br price weight price_* e*
	
	*-回归平方和 & 残差平方和
	regress price weight
	
	sum price
	gen u2=price_fit - r(mean)   
	egen ssr=sum(u2*u2)   //回归平方和=预测Y值与实际Y均值的差值的平方和
	format ssr %20.0f
	
	egen sse=sum(e*e)
	format sse %20.0f
				*-自由度df：以样本来估计总体时，样本中独立或能自由变化的个数
	* Model df：1表示回归模型有1个参数
	* Residual df：总自由度(n-1) - 模型自由度
	
	*-拟合优度 R-square (R2) & Adj R-squared 
	shellout "R1_OLS_venn.pptx"  //韦恩图, Venn
	
	corr price price_fit
	dis r(rho)^2
	
	/*	维基百科对Coefficient of determination（也就是R方）的解释：
	" In linear least squares multiple regression with an estimated intercept term, R^2  equals the square of the Pearson correlation coefficient between the observed Y and modeled (predicted)  Y_hat data values of the dependent variable." 
	
	在带有截距项的线性最小二乘多元回归中，R^2 等于实测值 Y 和拟合值 Y_hat 的Pearson相关系数的平方。     */
	
	regress price weight
    regress price weight trunk  //R2会受到变量个数的影响
	
	*R2_adj=1-((1-R2)(n-1)/(n-p)), 其中n为样本量，p为变量数
	
	*-OLS 的估计系数是一个随机变量, SE 衡量了其不确定程度;
	
	*- t = b/se
	regress price weight
	dis "t-value = " %4.2f _b[weight]/_se[weight]  // t = b/se

	
*-4.3.4- 常数项的含义

	use "B1_consume.dta", clear
	regress consume income
	
	*-Comments:
	* 1 常数项 = E(Y|x=0), 
	*   本例中，收入为零时的基本消费支出为 51.9 元.
	* 2 回归中遗漏常数项，相当于人为附加约束条件 constant=0
	*   换言之，相当于假设 E(Y|x=0)=0.
	  
	*-常数项的含义依赖于 x 的取值
	sysuse "nlsw88.dta", clear
	label var hours "hours"
	reg wage hours         // Model A
	est store A
	*-Note: 此时的常数项没有意义，因为 min(hours)=1
	*-      我们研究的是经济学问题，而不是纯粹的统计学问题
		
	sum hours
	gen hours_demean = hours-r(mean)   // (x-x_mean)
	label var hours_demean "hours-[Mean]"
	reg wage hours_demean  // Model B
	est store B
		
	gen hours_40 = hours-40            // (x-40)
	label var hours_40 "hours-40"
	reg wage hours_40      // Model C
	est store C
		
	esttab A B C, mtitle(A B C) nogap s(r2 N) label

*-4.3.5- 系数的含义

	*-在控制其他变量不变的情况下, x 对 y 的边际影响
	
	*-一个常见的 Puzzle
	sysuse "auto.dta", clear
	cap eststo clear
	eststo: reg price length
	eststo: reg price length mpg weight
	esttab, nogap

	shellout "R1_OLS_venn.pptx" 

*-4.3.6- 稳健型标准误

*-4.3.6.1	异方差稳健型标准误
	sysuse "auto.dta", clear
	reg price weight, robust // White(1980) 异方差稳健型标准误

	*-	例子
    sysuse "nlsw88.dta", clear 
    global x "hours ttl_exp tenure south collgrad married"
	
	*-去除离群值；以保证后续分析的样本数一致
	qui reg wage $x i.industry i.occupation
	keep if e(sample)
	
	*-Homo v.s. Het
	reg wage $x            // 干扰项同方差
	est store m1
	reg wage $x, robust    // 干扰项异方差
	est store m2
	esttab m1 m2, mtitle(homo robust) b(%4.3f) se(%6.4f) brackets ///
	star(* 0.1 ** 0.05 *** 0.01) nogap
	   
	*-Note：这是 90% 以上的文献都采用的方法；
	*       后续复杂模型的稳健型标准误也基本上以 White(1980) 为基础  

*-4.3.6.2	聚类调整后的标准误
	reg wage $x, vce(cluster industry)
	est store m3
	esttab m1 m2 m3
	*-Note: 
	* (1)   同行业内的个体的干扰项之间彼此相关;
	* (2) 不同行业之间的个体的干扰项之间彼此不相关；
	* (3) 系数估计值仍然采用OLS估计值，因为它是无偏的
	
*-4.3.6.3	自抽样法(Bootstrap)稳健型标准误
    use "B1_consume.dta", clear
    bootstrap, reps(300): reg consume income
	
	*-等价于  
	help vce_option
	reg consume income, vce(bs,reps(300))
	
    * 基本思想和步骤：
	*
    * 假设样本是母体中随机抽取的；
    * 通过反复从样本中随机抽取经验样本来模拟母体的分布；
	*
	* 1. 估计原始模型, 得到 x 的估计系数 b_x;
	*
    * 2. 从样本中可重复地(有放回地)抽取 N 个观察值，
	*    得到经验样本 S1(Empirical Sample)，对经验样本 S1 执行 OLS 估计，
	*    得到 x 的系数估计值 b_1；
    *
	* 3. 将第2步重复进行 k 次(如 k=500) 
	*    (相当于得到 k 组经验样本，S1, S2, ... Sk,
	*    每一组经验样本都可以视为对母体的一次随机抽样), 
	*    记录针对于每个经验样本得到的 x 变量的 OLS 估计值, 即
	*    b_j = {b1, b2, ..., b500}；
    *
	* 4. 统计这 500个估计值的标准差 sd(b_j) = sd{b1, b2, ..., b500}，
	*    将其视为实际估计值 b_x 的标准误, 即 sd(b_j) = se(b_x);
    *
	* 5. 计算 t 值: t = b_x/se(b_x), 以及相应的 p 值
	
	*-Note: 抽样次数与样本数和抽样目的有关;
	*       通常，抽样次数 K = 500 或 1000
	

*	----------------
*	更多知识
*	连享会：回归分析专题
	view browse "https://www.lianxh.cn/blogs/32.html"		




*			  ==========================================
*						     Stata公开课 
*  
*                           软件及计量基础   
*             ==========================================

*                    	 主讲人：候丹丹
*                   	 主办方：连享会（www.lianxh.cn）
*           课程主页：https://gitee.com/lianxh/stataopen
*                     :: 课件下载，答疑等 ::

*			       ============================
*					 	第四讲 回归分析（一）   
*                  ============================
*					    -4.4- 结果的呈现

*-4.4.1- logout

	ssc install logout,replace
	
	sysuse "nlsw88", clear
	global y "wage"        
	global x "hours tenure married collgrad" 
	
	*-去除缺漏值
    qui reg $y $x i.race i.union i.industry i.occupation
	keep if e(sample)
	
	*-输出描述性统计分析表
    logout, save("$Out\Table1_sum01") excel replace:  ///
            tabstat $y $x, column(stats) format(%6.2f) ///
            stats(mean sd min p50 max) 	
			
	*-输出相关系数矩阵
	logout, save("$Out\Table2_corr") excel replace: ///
            pwcorr_a $y $x
			
*-4.4.2-	esttab
	reg $y $x
	est store m1
	eststo m2: reg $y $x i.race
	eststo m3: reg $y $x i.ind
	local s "using $Out\myout.csv" //执行时，选中此行可以把结果输出到Excel
    esttab m1 m2 m3 `s', nogap replace  ///
	        s(N r2_a )  ///
	        sfmt(%6.0f %4.3f) 
	
*-4.4.2-	sum2docx
	ssc install sum2docx    //stata15版本以上可用
	
	*-输出描述性统计分析表
	sum2docx $y $x using sum2docx.docx, ///
		replace stats(N mean(%9.2f) sd(%9.2f) median p25 p75 min max) ///
		title("表1 描述性统计表")  ///
		note("this is auto's Summary Statistics")
	
	*-输出相关系数矩阵
	corr2docx $y $x using corr2docx.docx, ///
		replace star(* 0.05) fmt(%4.2f) ///.
		title("表2 相关系数矩阵")
		
	*-输出组间均值差异检验
	putdocx clear
	t2docx $y $x ///
        using t2docx.docx,   ///
        replace by(union)   title("表3 t检验")
	
	*-输出回归结果
	reg $y $x
	est store m1
	reg $y $x i.race
	est store m2
	reg2docx m1 m2 using result.docx,         ///
		ar2(%9.2f) b(%9.3f) t(%7.2f) r2(%9.3f) ///
        title("表4 回归结果") replace  
	
*-4.4.3.-	小技巧:回归结果中不报告虚拟变量的系数
	 
	sysuse "nlsw88", clear
	global y "wage"        
	global x "hours tenure married collgrad" 
	
	reg $y $x, r
	est store m1
	reg $y $x i.industry, r
	est store m2
	reg $y $x i.industry i.occupation, r
	est store m3
	
	local s  "using Table03.rtf"  // 输出到word文档
	local m  "m1 m2 m3"        // 模型名称
	local mt "model1 model2 model3" //模型标题
	esttab `m' `s', mtitle(`mt') b(%6.3f) nogap compress  ///
		star(* 0.1 ** 0.05 *** 0.01)  ///
		ar2 scalar(N) replace  ///
		indicate("Industry=*.industry" "Occupation=*.occupation" )
		





*	----------------
*	更多知识
*	Stata：回归结果中不报告行业虚拟变量的系数
	view browse "https://zhuanlan.zhihu.com/p/34940809"
*	知乎
	view browse  "https://www.zhihu.com/question/32021302"
*	33讲
	view browse    	"https://space.bilibili.com/546535876/channel/seriesdetail?sid=684350"	











